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We present a numerical study of an idealized magnetohydrodynamic (MHD) configuration con- 
sisting of a planar wake flow embedded into a three-dimensional (3D) sheared magnetic field. Our 
simulations investigate the possibility for in-situ development of large-scale compressive disturbances 
at cospatial current sheet - velocity shear regions in the heliosphere. Using a linear MHD solver, we 
first systematically chart the destabilized wavenumbers, corresponding growth rates, and physical 
parameter ranges for dominant 3D sinuous-type instabilities in an equilibrium wake-current sheet 
system. Wakes bounded by sufficiently supersonic (Mach number Ms > 2.6) flow streams are found 
to support dominant fully 3D sinuous instabilities when the plasma beta is of order unity. Fully 
nonlinear, compressible 2.5D and 3D MHD simulations show the self-consistent formation of shock 
fronts of fast magnetosonic type. They carry density perturbations far away from the wake's cen- 
ter. Shock formation conditions are identified in sonic and Alfvenic Mach number parameter space. 
Depending on the wake velocity contrast and magnetic field magnitude, as well as on the initial 
perturbation, the emerging shock patterns can be plane-parallel as well as fully three-dimensionally 
structured. Similar large-scale transients could therefore originate at distances far above coronal 
helmet streamers or at the location of the ecliptic current sheet. 

PACS numbers: 52.30.Cv, 52.65.Kj, 96.50.Ci, 96.50.Fm 



I. INTRODUCTION 

The magnetic and plasma flow configuration in a solar 
coronal streamer belt beyond its helmet cusp is character- 
ized by a region of cospatial magnetic shear and a cross- 
stream velocity variation reminiscent of a 'wake' flow. 
Suitably simplified, a local study of streamer belt dynam- 
ics may then consider a 'boxed' model of a fluid wake - 
current sheet configuration as sketched in Fig.Q] A large 
number of previous studies have identified the various lin- 
ear MHD instabilities supported by such an equilibrium 
configuration, with purely incompressible studies as 
well as extensions to the compressible regime 0, 0,M 1^ ■ 
Their applications cover a range of physical problems, as 
they are used to model dynamics at the heliospheric cur- 
rent sheet 0, mechanisms for variable slow solar wind 
formation and acceleration Jj, and plasmoid formation 
and evolution in the solar streamer belt 0, 0, • 

Similarly idealized configurations equally containing 
velocity shear layers co-spatial with an electric current 
sheet or neutral sheet have been used to investigate the 
linear and nonlinear stability properties of planar and 
non-planar current - vortex sheets 8, 9^(I11. 12i. [T^[l^ 
and planar magnetized jets llJillE E3 13 ■ These con- 
figurations were investigated in the context of solar flares 
and surges, among other applications. The possibility 
of having both shear-flow related instabilities, and ideal 
and resistive magnetic instabilities, allow the plasma to 
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transit to turbulent regimes via various routes in its fur- 
ther nonlinear evolution. Energy can be tapped from 
the background mean flow, with localized reconnection 
events releasing magnetic energy. 

In this paper, we revisit the stability properties for a 
compressible wake - current sheet model of streamer belt 
dynamics, and extend earlier two-dimensional nonlinear 
simulations 0, Q ^3 to three-dimensional MHD compu- 
tations. The necessity to model the compressive aspects 
when applying the 'box-model' from Fig.^to encompass 
both the slow and fast solar wind regions is obvious since 
the wind becomes supersonic sufficiently close to the Sun 
(at roughly 5Rq, as follows from Parker's model [T^). 
Compressibility can give rise to new effects as compared 
to incompressible studies, such as nonlinear steepening 
of ideal sinuous instabilities under supersonic flow con- 
ditions to form fast magnetosonic shock fronts 0, 
and the domination of fully three-dimensionally struc- 
tured unstable modes over two-dimensional ones 0, • 
In the incompressible case, the validity of Squires theo- 
rem Q precluded the latter possibility and allowed one 
to concentrate on 2D dynamics. It is the main result 
of this paper to show that (1) strongly supersonic mag- 
netized wakes, characterized by a plasma beta of order 
unity, also have dominant 3D sinuous instabilities which 
(2) can lead to fully 3D structured fast magnetosonic 
shock fronts in their nonlinear evolution. 

Earlier, formation and ampliflcation of fast magne- 
tosonic shocks was shown for supersonic super-Alfvenic 
wakes by means of numerical simulations 0- These 
shocks were also observed in studies of 2D current-vortex 
sheets in the transonic regime llj. We supplement these 
earlier findings by a complete parametric study of ini- 
tial state parameters which lead to shock formation. 
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Our findings that dominant 3D sinuous instabifities are 
possible correct earfier claims to the contrary where 
wavenumber space was not suflticiently explored. At the 
same time this calls for a study of the nonlinear develop- 
ment of a sinuous instability in three dimensions, which 
we perform as well. 

It was pointed out by Wang et al. on the basis 
of a linear stability analysis of a similar wake - current 
sheet configuration, that instabilities may develop in-situ 
at the heliospheric current sheet at radial distances of 
order 0.5 - 1.5 AU. They emphasized the importance 
of streaming sausage modes leading to traveling mag- 
netic islands (clouds) and plasmoids. Our results focus 
on the streaming kink mode, and our nonlinear simula- 
tions demonstrate that it can lead to shock-dominated 
large-scale transients. This re-emphasizes the possibil- 
ity that not all satellite solar wind observations of such 
events need to be directly mapped back to a specific solar 
coronal event (eruptive prominence or CME). 

This paper is organized as follows. Section ITU presents 
the system of MHD equations describing the wake- 
current sheet evolution and the initial equilibrium con- 
figuration studied. Section IIIII describes the numerical 
tools used for the linear stability analysis and for the 
nonlinear evolution studies. In Section lWI we present the 
linear stability analysis of the wake-current sheet config- 
uration and identify the previously unexplored domain 
of the dominating three-dimensional sinuous instability. 
In Section the nonlinear MHD simulations with the 
Versatile Advection Code J_9] for different values of wake 
velocity contrast and magnetic field magnitude are dis- 
cussed. Formation of fast magnetosonic shocks is demon- 
strated and their properties and generation conditions 
are described. Conclusions and a discussion of the rel- 
evance to in-situ streamer belt dynamics are placed in 
Section Ell 



II. MHD EQUATIONS AND EQUILIBRIUM 
MODEL 

The system under consideration is described by the set 
of one-fluid compressible resistive MHD equations of the 
form: 

dtp + V- (pv) - 0, 
dt (pv) + V • (vpv - BB) + Vptot = 0, 
^te -h V • (ve - BB ■ V + vptot) = V • (B x 77J) , 
5tB -I- V • (vB - Bv) = -V X (ryJ) , 
J = V X B, 

p = (7 - 1) (e - pvy2 - BV2) , Ptot =P + BV2, 

where p is the mass density, v is the plasma velocity, B 
is the magnetic field, and e stands for total energy den- 
sity. Other quantities are J - electric current, thermal 
pressure p and the constants 7 = 5/3- adiabatic index, 
and 7] - resistivity coefficient. The resistive source terms 



are taken along to accurately simulate the small-scale re- 
connection events which only occur far into the nonlinear 
regime in the cases studied below. 

Guided by the box model of Fig. ^ the magnetized 
wake fiow co-spatial with the current sheet can be mod- 
eled as a force-free equilibrium configuration: 

V,^l- cosh-iy, Vy^Q, K = 
= ACHanh-, B,, ^Q, B^^ ilC^cosh"^- 

p=1.0, p=(7M2)-i, (1) 

where Ms and Ma are the sonic and Alfvenic Mach num- 
bers for the fast flow streams, respectively, and w de- 
scribes the thickness of the current sheet relative to the 
width of the wake flow. The latter, together with the 
density and velocity of the fast flow streams, has been 
used to define our unit system. The plasma /? is then 
found from /3 = (2M^)/(7Aff ). If not mentioned oth- 
erwise, we set w = 1. The same force- free equilibrium 
was analyzed in 2], while a similar, but non force-free 
magnetic configuration, was used in 0]. 



III. NUMERICAL TOOLS 

Linear stability analysis of the wake-current sheet 
configuration was performed numerically with the 
LEDAFLOW code This code can compute the 

complete MHD spectrum of all waves and instabilities 
that are eigenfrequencies of one-dimensionally varying 
equilibria containing background plasma fiows. The lin- 
earized MHD equations are discretized using finite ele- 
ments (quadratic and cubic Hermite polynomials) in the 
direction of inhomogeneity and a Fourier representation 
of the perturbations in the invariant directions. Assum- 
ing a temporal variation as exp(At), application of the 
Galerkin procedure results in a complex non-Hermitian 
eigenvalue problem which is then solved by a complex 
QR-decomposition-based solver. For each isolated eigen- 
value A, the set of corresponding eigenfunctions can be 
found with an inverse vector iteration. In the QR mode, 
the whole spectrum for a given set of Fourier mode num- 
bers is calculated, but the amount of grid points in the 
finite element representation is usually limited to 100 - 
120 points. In inverse iteration mode, the spatial resolu- 
tion is much higher to be able to get converged eigenfre- 
quencies with strong gradients in their eigenfunctions (up 
to 5000 points). The boundary conditions implemented 
in LEDAFLOW assume perfectly conducting walls at 
y-boundaries, which is, strictly speaking, not applica- 
ble to the streamer belt configuration. To prevent any 
influence of these conducting walls on the solution of 
the eigenvalue problem, one should place them at suf- 
ficiently large distances in the perpendicular y-direction, 
and check whether the eigenfunctions are localized away 
from the computational domain boundaries and vanish 
close to the walls. 
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Nonlinear modeling of the wake-current sheet dynam- 
ics was performed by the general finite- volume-based Ver- 
satile Advection Code (VAC) [ll|2l|. For all runs we 
used the full set of compressible resistive MHD equa- 
tions with a constant resistivity 77 — 10~* value. The 
one-step total variation diminishing (TVD) method |23 
was employed, which is a second-order accurate shock- 
capturing scheme. It uses an approximate Roe-type 
Riemann solver with the minmod (for 3D simulations) or 
Woodward (for 2.5D simulations) limiting applied to the 
characteristic waves. 

The rectangular computational domain has periodic 
boundaries in the streamwise x-direction, and open 
boundaries in the cross-stream direction y, such that 
perturbations reaching these boundaries are allowed to 
propagate freely outside the computational box. For 3D 
runs, the spanwise direction z is periodic as well. We 
used non-equidistant rectangular grids with a symmetric 
(around y = 0) accumulation near the core of the wake. 
Grid resolutions varied from 100 x 200 up to 300 x 600 in 
2.5D and from 60 x 90 x 60 up to 100 x 150 x 100 in 3D 
cases. The grid lengths in the streamwise and spanwise 
directions are found from — 27: /kx and — 
respectively, with wavenumbers kx, kz corresponding to 
the most unstable mode identified from the linear sta- 
bility analysis. For the cross-stream boundary, the open 
boundaries are placed at Ly = ±25-50 depending on the 
localization of the unstable eigenmode. 2.5D runs were 
performed on an SGI Octane workstation, while 3D sim- 
ulations were executed on typically 32 processors of the 
SGI Origin 3800 at the Dutch Supercomputing Centre, 
using openMP for parallelization purposes. 



IV. LINEAR STABILITY ANALYSIS 

Depending on the parameters of the one-dimensionally 
varying (but fully three-dimensionally structured) equi- 
librium configuration given by Eqs. (Q, as well as on the 
value of the resistivity and the lenghtscales of the im- 
posed perturbation, up to three different mode types may 
be destabilized in the wake-current sheet configuration. 
These are known as the ideal sinuous (or streaming kink) 
mode, and an ideal and a resistive varicose (sausage- type) 
mode 0, 01 ■ In a 2D case these types are easily distin- 
guishable by the symmetry in their cross-stream velocity 
perturbations, which is even for the sinuous mode and 
odd for both varicose modes. In a 3D case these symme- 
try properties are broken. Here we shall mainly focus at 
the instability of an ideal sinuous type, already present 
in unmagnetized hydrodynamic wake configurations. In 
the range of Mach and Alfven Mach numbers investi- 
gated, its growth rate is always larger than that of the 
varicose modes, and thus the sinuous mode is expected to 
dominate the evolution of the system. This mode is also 
known to demonstrate interesting nonlinear evolution: in 
2.5D MHD simulations, sinuous mode instabilities led to 
the self-consistent formation of fast magnetosonic shocks 



in the nonlinear stage 

In Fig. 121 a typical spectrum of a wake-current sheet 
system obtained with LEDAFLOW is presented. The 
parameters used are for a supersonic = 3, super- 
Alfvenic Ma = 5 wake at fixed streamwise wavenumber 
kx — 0.35 and spanwise wavenumber k^ — 0.15. The fig- 
ure shows individual complex eigenfrequencies A through- 
out the complex plane, but for the value of the resistivity 
rj = 0.01, the ideal and resistive varicose instabilities are 
fully suppressed, and only the ideal sinuous mode re- 
mains unstable in the system (marked by the arrow in 
Fig. 121) . Since the sinuous instability is ideal in nature, 
its growth rate does not depend on the value of the resis- 
tivity chosen. The eigenfunctions, namely the perturbed 
plasma density Sp and cross-stream velocity 6vy , for the 
unstable sinuous mode are shown in Fig. [31 The insta- 
bility growth rate is 7 « 0.0402 and the phase velocity 
Vph ~ 0.55. It is seen that 6vy is nearly symmetric (even) 
with respect to y = 0, while the density perturbation has 
a more complicated structure but is in essence odd. All 
linear stability results following are performed for a fixed 
value of resistivity as used in the nonlinear simulations, 
namely t] = 10^'*. 

The maximal sinuous mode growth rate is achieved for 
an incompressible (limit AIs — > 0), purely hydrodynami- 
cal (B=0) case. In incompressible magnetized wakes, the 
dominant instability was found to be two-dimensional 
having zero spanwise wavenumber, in agreement with 
the Squires theorem 0. When plasma compressibility 
is taken into account, as done in this paper, the situa- 
tion can change drastically. In [3], it was demonstrated 
that at high Mach numbers, an increase of the spanwise 
wavenumber can make the ideal compressible varicose 
mode more unstable. The sinuous mode growth rate, 
however, was believed to decrease with increasing k^. 

In Fig. ^ the numerically obtained growth rates for 
Ms = 3, Ma = 5 wakes of the ideal sinuous instability 
are plotted versus streamwise wavenumber k^ (left panel 
A), and spanwise wavenumber k^ (right panel B). The 
top curve in panel (A) for fc^ = is identical to the 
corresponding dispersion curve for Mach 3 wakes shown 
in 2], their Fig. 1(b). It is clearly seen from panel (A) 
that the region of destabilized kx is widest in the 2D case 
(kz — 0). Increasing the spanwise wavenumber kz leads 
to a shrinking of the instability domain from both sides, 
while for increasing kz, the most unstable streamwise 
wavenumber also increases. This trend continues up to 
some upper value of kz > 0.6 where the sinuous mode be- 
comes stable. Panel (B) demonstrates that the 2D land- 
scape of growth rate versus streamwise kx and spanwise 
kz wavenumber, has a secondary 'ridge' which was over- 
looked in 0]. There, at a fixed wavenumber kx = 0.55, an 
increase of the spanwise wavenumber led to a decreasing 
growth rate for the sinuous mode (their Fig. 3(b)). Our 
panel (B) shows that their study overlooked the rather 
interesting case of larger streamwise wavelength (smaller 
kx), where three-dimensional local extremum of growth 
rate was found for kx < 0.25. We can further observe 
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that with an increase of the streamwise wavenumber 
from 0.075 up to values 0.35, both the growth rate and 
the instabihty existence domain (in kz) increase, reaching 
their maximum values at k^ ~ 0.35, which corresponds 
to maximal 2D growth. Further increasing k^ leads to 
a renewed drop in growth rate, however the instability 
domain does not change significantly. The latter can be 
seen when comparing our = 0.35 curve from panel (B) 
with the kx = 0.55 curve for Mach 3 wake flow in 0, their 
Fig. 3(b). Below some kx < 0.075 the sinuous mode is 
fully stabilized again. 

The importance of this secondary ridge in the growth 
rate landscape is marginal for the Mg = 3, Ma = 5 
flow conditions of Fig. ^ Indeed, the absolute maxi- 
mal growth rate corresponds to a purely 2D {k^ — 0) 
mode with kx ~ 0.35. However, for different combina- 
tions of Mach and Alfven Mach number, the most un- 
stable mode can be fully 3D in nature, a fact completely 
missed by earlier studies of this wake system. Variations 
of ideal sinuous mode maximal growth rate versus span- 
wise wavenumber for different values of Mg and Ala are 
plotted in Fig. O and the possibility of a 3D dominant 
instability is clearly manifested. Note that in this fig- 
ure, the streamwise wavenumber kx varies from curve to 
curve, and has a value corresponding to the overall max- 
imum in the growth rate landscape. 

For a given Mach number, increasing the magnetic 
field strength (decrease of Alfven Mach number Ma) acts 
to suppress the sinuous instability. This is also seen in 
Fig- El panel (A). The range of destabihzed kx decreases, 
and the streamwise wavelength of the fastest growing 
mode decreases. When Ma becomes smaller than some 
critical value M^'^, the system is stable. This critical 
value Ala^ was found to be approximately equal to 2.5, 
and this value is rather universal in a wide range of stud- 
ied cases (0 < Mg < 10). This is consistent with ^j, their 
Fig. 2. Note that Wang et al. Q found a lower cut-off 
value (namely M^*" ?s 1.2) which must be due to the dif- 
ferent equilibrium configuration used. Their equilibrium 
magnetic field strength is weaker at the center of the flow 
wake, while our configuration has a constant field magni- 
tude throughout. Therefore, the stabilizing influence of 
magnetic tension translates into a lower cut-off value for 

M^r 

Generally, with an increase of the sonic Mach number 
Ms , the overall maximal growth rate of the sinuous insta- 
bility goes down, as well as the diapason of destabilized 
wavenumbers. For an Alfven Mach number Ma = 3.5, 
variations from Mach 2 till Mach 15 flows on the sinuous 
growth rates are shown in panel (B) from Fig. El Corre- 
sponding dependences TZe [X{kx)] can be found elsewhere, 
e.g. see Fig. 1(b) in 0. 

The most important finding from Fig. El is that the 
most unstable sinuous mode may become oblique to the 
shear layer in a wake - current sheet system for suffi- 
ciently supersonic flows. Untill now, this was known for 
the compressible plane current-vortex sheet [l^ as well 
as for varicose- type modes in compressible wake - current 



sheet configurations We found that for a given sonic 
Mach number larger than Mg « 2.6, there exists a range 
in Ma where the maximum growth occurs for a mode with 
nonzero spanwise wavenumber. That diapason of Ma for 
dominant three-dimensional instability is bounded from 
below and from above, and most wide near Mg ~ 2.6 
where it is 2.9 < Ma < 4. Its size decreases slowly with 
increasing Mg, so that even for highly supersonic wake 
flows the 3D instability domain is non-negligible, e.g. for 
Mg = 10 it is found in the range 2.65 < Ma < 3.025. 

In the next sections we will study the nonlinear evo- 
lution of the wake-current sheet system in 2D and 3D 
using the information about the most unstable modes as 
calculated by LEDAFLOW for constructing the initial 
conditions. 



V. NONLINEAR EVOLUTION AND SHOCK 
FORMATION 

A. 2.5D MHD parameter study of shock formation 

Since we are mainly interested in the nonlinear evolu- 
tion of sinuous- type modes which have the largest growth 
rates and therefore dominate the system, the initial equi- 
librium configuration given by Eqs. was perturbed 
by a symmetric cross-stream velocity perturbation of the 
form: 

Svy ~ £ sin (kxx) exp(— y^), 

with typical values of e = 10^'^..10~^. A comparative 
analysis of the nonlinear wake response to such per- 
turbations in 2.5D MHD, for varying sonic and Alfven 
Mach numbers, was presented in j2i- Here, we will focus 
on previously unaddressed issues, namely the interplay 
between simultaneously excited modes having different 
lenghtscales and growth rates, and the properties and 
formation conditions for shock waves, which are observed 
in the simulations at the instability saturation stage. 



Mode-mode interplay 

In Fig. El an example of a 2.5D MHD simulation 
is presented. By definition, all 2.5D simulations in- 
clude only kz = modes. Initially, two sinuous modes 
with streamwise wavenumbers kxi = 0.35 (mode 1) and 
kx2 = 0.35/3 = 0.11667 (mode 2) were excited on the grid 
with amplitudes ei — 2.5 ■ 10~^ and 62 = 10"^, respec- 
tively. For the simulated Mg — 3, Ma — 5 {(3 — 3.33) 
wake, mode 1 corresponds to the dominant 2D fastest 
growing mode in the system. Note that we perturb the 
slower growing mode at a higher amplitude to force mode 
interaction. Fig. El shows snapshots of density and cross- 
stream velocity, while for the same time frames. Fig. [3 
presents the spatial spectrum of the cross-stream veloc- 
ity. This spatial Fourier transform performed for the 
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sequence of snapshots allows to extract information on 
the linear growth rate for each mode and to confirm 
LEDAFLOW predictions. The linear results of the previ- 
ous section were found to be very accurate, the difference 
in growth rates deduced from VAC simulations and those 
obtained by LEDAFLOW is usually less than 1%. 

At early times, mode 2 is seen to grow with the growth 
rate « 0.018. However, at times around T = 200, as it 
is clearly seen in the spectral evolution, mode 1 obtains 
a comparable amplitude. Having the larger growth rate, 
mode 1 develops faster, and at times around T = 250 its 
amplitude becomes larger than that of mode 2. Later on 
in the nonlinear evolution the short-wavelength mode 1 
remains the most dominant one, but a cascade to smaller 
wavelength dynamics is obviously occuring. Note that 
initially (times T < 200) this only involves direct over- 
tones of the two excited modes. 

Similar observations can be made as follows. In the 
linear stage of the instability, the volume average cross- 
stream kinetic energy in the system plotted in Fig.|Hl and 
defined as 

Ehn = {2L,LyL,)-^ / / / pvldxdydz, 

Jo J-0.5Ly Jo 

(2) 

grows exponentially in accordance with the slower grow- 
ing mode 2. Again, its slope changes near T = 200 in- 
dicating that the faster growing mode 1 comes into play. 
Eventually, transverse kinetic energy growth is halted 
and this saturation point typically coincides with the 
time where the largest density contrast is achieved. The 
total kinetic energy of the system decreases, while the 
magnetic energy increases indicative of the occuring en- 
ergy conversion processes. 



Tresholds for shock formation 

Depending on the system parameters, exactly at the 
point of maximal density contrast, shock waves may be 
formed. They are clearly detected in the density snap- 
shot at time T — 300 from Fig. |S1 Note that these shocks 
extend far out into the fast stream regions, up to 25 
wake widths in the cross-stream direction. Fig. |51 shows 
a cross-cut of the density variation at this time, taken at 
about 10 wake widths above the wake-current sheet. The 
shocks carry strong density contrasts, up to a factor of 
2. Note that due to the mode-mode interaction, succes- 
sive shock fronts can have different strengths, and their 
interseparation no longer corresponds to the wavelength 
of the linearly most unstable mode. This cross-cut can 
loosely be interpreted as a time-evolution of an in-situ 
satellite measurement, since the shock fronts are typi- 
cally advected over the periodic side boundaries for more 
than a full advection cycle. 

To identify the type of the shock front, we used 
Rankine-Hugoniot relations to determine the shock 
speed. This is most easily found from mass continu- 



ity, where the velocity of the shocks is found as Vsh = 
{{pvi^t^"''" - {pv^ff''')/{p*"'f°^^ - The shock 

type then follows from determining Alfvcn, slow- and fast 
magnetosonic Mach number transitions in the coordinate 
system which moves with the shock velocity. Shock waves 
formed at the nonlinear stage of the sinuous instability 
of a supersonic magnetized wake were found to be fast 
magnetosonic in nature. 

The steepening of the sinuous wave fronts into shock 
waves is possible only above a threshold value of sonic 
Mach number. That threshold depends on the ampli- 
tude of the magnetic field (hence on Ma). In a pure hy- 
drodynamical case {Ma — oo), sonic shocks are formed 
for supersonic wake flows with Mg > 2.10. Shocks ap- 
peared for all studied hydrodynamic cases in the range 
2.10 < Ms < 10.0 and no upper hmit of Mg for their 
formation was found. This hydro critical value is recov- 
ered in weak {Ma > 5) magnetized wakes. However, 
with an increase of the magnetic field magnitude, the 
minimum value of the sonic Mach number needed 

for shock formation, is increased. This is presented in 
Fig. ^1 where threshold sonic Mach numbers are plot- 
ted versus the Alfven Mach numbers. Formation of the 
shock fronts was shown to be possible for any Ma where 
the sinuous mode is still unstable {Ala > 2.5). Hence, 
for any magnetic field strength, where the sinuous insta- 
bility is not suppressed linearly, we were able to iden- 
tify a shock-treshold value of Mach number M^™ . For 
all Ms > Mf'"'{Ma) shock fronts form self-consistently 
through wave front steepening. 

Most of the nonlinear runs to determine the parame- 
ter ranges for shock formation were performed for sys- 
tems of length Lx — 27r/fc™, with fc™ being the stream- 
wise wavenumber of the fastest growing mode. It turns 
out that shocks can appear having different scales in the 
x-direction, since there exists a diapason of streamwise 
wavenumbers kx, located around /c™, for which shocks 
develop. At scales sufficiently far from fc™, formation 
of shock waves was not observed. For a given Alfven 
Mach number, with an increase of the wake speed con- 
trast {Ms), the angle between the shock front and the 
streamwise direction decreases. Hence the density per- 
turbations become more and more cuddled up to the core 
of the wake, and thus not always reach far cross-stream 
regions. 

Finally, when applied to the coronal streamer belt, the 
width of the current sheet - w in Eqs. iQJ - is typically 
much smaller than that of the fluid wake 3. With a 
decrease of the relative width of the current sheet, the 
number of grid points needed to achieve acceptable reso- 
lutions grows rapidly, and simulations become more and 
more time- and memory consuming. Nevertheless, for 
selected cases we did vary the current sheet width, and 
these experiments did not show any qualitative differ- 
ence in the nonlinear stage of evolution, indicating that 
the model of equal-sized wake and current sheet used 
throughout this paper is warranted. Quantitatively, as 
can also be obtained from a linear LEDAFLOW analy- 
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sis, the growth rate and streamwise wavenumber of the 
most unstable sinuous mode decreases when the current 
sheet width goes down. 

B. Case studies in 3D MHD 

The 2.5D simulations in the previous section excluded 
the emergence of fully 3D structured shock fronts. How- 
ever, our linear stability study demonstrated the exis- 
tence of dominant 3D unstable modes. Here we discuss 
three particular three-dimensional runs performed for dif- 
ferent parameters of the equilibrium configuration and 
for different initial excitations. We investigate two dis- 
tinct cases: (1) a wake with a dominant purely 2D mode, 
where we look into the possibility of 3D shock structur- 
ing due to mode-mode interactions; and (2) a wake pre- 
viously identified to have a dominant 3D linear mode. 



Mode-mode interaction 

In Fig. 111! plasma density isosurfaces are plotted for 
different stages in the sinuous mode instability develop- 
ment. In the first column of Fig. values of sonic 
and Alfven Mach numbers were chosen to be Mg = 3, 
Ma = 5, as in the 2.5D run from Fig.|Hl The wake was 
perturbed with two modes: a pure 2D mode - the most 
dominant instability - having = 0.35, = 0, and 
a 3D mode k^ = 0.35/3 = 0.11667, = 0.125. The 
amplitude of the 3D mode (10~^) was much larger than 
that of the 2D mode (10~^). At first, the 3D sinuous 
instability develops in accordance with its growth rate 
w 0.018, and the emerging density structure is clearly 
three-dimensional. However at times around T « 175 
(middle snaphot), the 2D mode comes into play. It leads 
to elongation of the density structures along the z-axis. 
It again causes a slope change in the Ej,^^ (t) dependence 
plotted in Fig.^] This continues untill the wake density 
structure tends to be mainly two-dimensional, with very 
slight modulations along the z-axis. In particular, the 
emerging shocks are nearly plane-parallel, as in the 2.5D 
simulations. 

If the same system is initially perturbed by the 3D 
mode alone, as is done in the middle column in Fig. 1111 
one can obtain fully 3D structured shock-dominated wake 
dynamics. The two-dimensional mode starts to grow 
from the noise level in this case, and obtains noticeable 
amplitude only at times when the 3D instability is close 
to saturation. The shock fronts have a truly 3D topol- 
ogy. Note that at times around T — 250 each high den- 
sity blob splits into parts whose scale corresponds to the 
scale of the fastly growing 2D k^ = 0.35, kz = mode. 
After saturation, a set of nearly two-dimensional (as seen 
from the third snapshot in the middle column in Fig. lll() 
structures form. The number of structures corresponds 
to the scale of the initially excited 3D mode, but their 
size in the streamwise direction is strongly influenced by 



the kx = 0.35, fc^ = mode. These two runs demon- 
strate that 3D shock structuring can even emerge from 
wakes with dominant 2D linear instabilities, but only un- 
der rather selective initial excitations. 



Dominant 3D instability 

The last column in Fig. 1111 shows the nonlinear evo- 
lution of the density in pure three-dimensional dynam- 
ics of the sinuous instability for the case Mg = 3, 
Ma — 3.5 {[3 — 1.6). Above, it was demonstrated that 
for these parameters, the most unstable sinuous mode 
is oblique to the wake's core and has the wavenumbers 
kx = 0.375, kz = 0.15 (see Fig. EJ. Note how at all 
stages in the evolution, the spatial density distribution 
remains fully 3D. As compared to the other cases shown, 
the maximal density contrast reached in the nonlinear 
stages is larger. For these parameters, the fast magne- 
tosonic shock fronts emerging in the nonlinear stage have 
a fairly complex 3D structure. The transverse kinetic en- 
ergy evolution is plotted in Fig. ^1 



VI. CONCLUSIONS 

In this paper we performed a numerical study of the 
linear properties and nonlinear evolution of wake ~ cur- 
rent sheet configurations by means of compressible resis- 
tive MHD simulations. We focused on three-dimensional 
effects for the ideal sinuous-type instability. In contrast 
to previous studies, we have shown that there exist ranges 
of Mach number Mg and Alfven Mach number Ma, where 
the most unstable sinuous mode is oblique to the shear 
layer. Such pure three-dimensional instabilities develop 
nonlinearly into fully 3D deformations of the wake sys- 
tem containing fast magnetosonic shock fronts. For dom- 
inant 3D instabilities, the Mach number should exceed a 
threshold value Mg > 2.6. With increase of Mg, the 
range of Ma corresponding to 3D instability decreases, 
but even for highly supersonic cases, a dominant 3D in- 
stability still exists in some range of Ma ■ At the nonlinear 
stage of this sinuous instability in the supersonic regime, 
self-consistent formation of fast magnetosonic shocks is 
observed. They carry density jumps (by factors of 2) 
in cross-stream fast flow regions far away from the wake 
center. As a result of mode-mode interaction, the shock 
strengths and interseparations can vary. Depending on 
the wake velocity contrast and magnetic field magnitude, 
as well as on the streamwise and spanwise scale of the ini- 
tial perturbation, these shocks may have plane-parallel as 
well as complicated three-dimensional structure. Shock 
fronts appear only for Mach numbers above some critical 
value, which depends on magnetic field strength, but is 
always larger than the hydro limit Mg = 2.1. 

Finally, let us estimate the heliospheric regions, where 
our box-model predicts the possibility for in-situ com- 
pressible transients. Using the Parker model ^Igj for the 
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interplanetary magnetic field one finds, in a manner orig- 
inally done by Wang et al. Alfven and sonic Macli 
number radial dependencies of the form: 



M,,(r) 



6v{r) 



2 



Mair) 



Sv{r) 



2r2 



1' 



where Sv{r) is the velocity contrast between fast and 
slow wind streams, and /3o and indicate the plasma 
beta and Alfven speed at 1 AU. For solar wind pa- 
rameters at the Earth's orbit given by n = 5 cm~'^, 
B = 5nT, Te = 2 • lO^K, Tp = 4 • lO^'K (taken from Q), 
one finds /3o ~ 1-6, while « 50 km/sec. Assuming 
that the fast-slow solar wind velocity contrast, as deter- 
mined by Ulysses measurements |25l | to be of the order 
dv « 350km/sec, does not change significantly with ra- 
dius beyond the various magnetosonic transitions, one 
can compute Ms{r) and Ma{r) dependencies. 

Our linear theory predicted that the sinuous mode is 
destabilized for Ma > 2.5 which then corresponds to dis- 
tances beyond 0.3 AU from the Sun. Relevant helio- 
spheric distances turn out to be 0.3 AU - 1.5 AU. Since 
we used the fluid wake scale length and the fast-slow ve- 
locity contrast to define dimensionless variables, we can 
verify whether the instability grows sufficiently fast to 



influence conditions at 1 AU. The flow wake width is ob- 
servationally determined to be of order 300.000km at 1 
AU, while at a heliospheric distance of order 0.5 AU it is 
estimated to be around 150.000km '24*1 . Our time unit 
will then correspond to approximately 7 minutes. At 0.5 
AU, the most unstable sinuous mode has a dimension- 
less growth rate of 0.0272 and corresponding streamwise 
wavenumber of fc^: = 0.325. Hence, the sinuous instability 
could fully develop within about 25 hours, in close cor- 
respondence with the transit time to 1 AU. The length- 
scale of these compressive perturbations are then of order 
2 • lO^^AU. This points to the distinct possibility of in- 
situ shock formation at wake-current sheet configurations 
within the heliosphere. 
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Computational domain 



FIG. 1: Model of the coronal streamer belt: a schematic view. 
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FIG. 2: MHD spectrum of a compressible wake-current sheet 
with Ms = 3, Ma = 5, fc^ = 0.35, fc^ = 0.15, 77 = 0.01. The 
unstable sinuous mode is marked. 



FIG. 3: Eigenfunctions corresponding to the unstable sinuous 
mode from Fig. |21 Real and imaginary parts of the cross- 
stream perturbed velocity and perturbed density are shown. 
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FIG. 4: Ideal sinuous mode growth rate versus stream- 
wise wavonumber kx and spanwise wavenumbcr kz- (A) 
At fixed spanwise wavenumber kz = {0.0, 0.25, 0.5, 0.6} top 
to bottom. (B) At fixed strcamwiso wavenumber kx = 
{0.35, 0.25, 0.15, 0.11667, 0.1, 0.075} top to bottom. 
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FIG. 5: Maximum growth rates of the ideal sinuous mode in a 
wake - current sheet system versus spanwise wavenumbcr for 
different Ms and Ma. Note that the strcamwise wavenumber 
kx varies from curve to curve. (A) At fixed Ms = 3, curves 
corresponding to Ma = {5.0,4.5,4.0,3.5,3.25,3.0,2.9,2.8} 
top to bottom, (B) At fixed Ma = 3.5, curves correspond- 
ing to Ms = {2.0, 2.5, 3.5, 4.0, 5.0, 15.0} top to bottom. 
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FIG. 6: 2.5D MHD simulations of a wako-currcnt sheet con- 
figuration with Ms = 3, Ma = 5, Tj = fO^*. Plasma density 
(left panels) and cross-stream velocity (right panels) at times 
T = 100, 175, 200, 250, 300, 500. Two sinuous modes were ex- 
cited and shocks of different strengths emerge at T = 300. 
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FIG. 7: Fourier spectrum of cross-stream velocity correspond- 
ing to the 2.5D run from Fig. |H| 



14 



FIG. 8: Vertical kinetic energy, total kinetic energy and mag- FIG. 10: Threshold value of sonic Mach number necessary for 
netic energy versus time for the 2.5D run from Fig.|^ shock formation versus Alfven Mach number. 
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FIG. 9: Plasma density distribution along the streamwise di- 
rection, far into the fast stream regions (y = 10.1) at T = 300, 
as obtained from Fig.|S| The fast shocks carry density jumps 
up to factors of two. 



3 


4 




3 


2 




3 







2 


8 




2 


6 


\ Shocks 


2 


4 




2 


2 


'No shocksV^^ ; 
HD limit^"^* 


2 








2 3 4 5 6 7 
M 



15 



FIG. 11: 3D MHD simulations of wake-current sheets. First 
column: Ms — 3, Ma = 5, two modes excited: {k^ = 
0.35, k-, = 0) and (fc^ = 0.35/3 = 0.11667, k^ = 0.125). 
Second column: as in first column, only the 3D mode ex- 
cited initially. Third column: Ms = 3, Ma = 3.5, ini- 
tiated with its dominant three-dimensional sinuous mode 
{kx = 0.375, fcz = 0.15). Density isosurfaces at consecutive 
times (top to bottom) are shown. 



FIG. 12: Evolution of vertical kinetic energies for all three 3D 
runs in Fig. 1111 



